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This  research  progressed  simultaneously  in  several  directions.  Details  of  each  area  of 
research  are  discussed  below. 

First,  we  developed  a  novel  neural  network  based  approach  to  imaging  ionospheric 
motion.  The  algorithm  employs  neural  networks  as  estimators  of  a  time  varying  electron 
density  distribution;  however  the  networks  are  not  used  to  estimate  the  electron 
densities  directly  but  rather  the  change  in  electron  densities  at  each  time  step. 

At  each  time  step,  the  method  consists  of  four  main  parts.  First,  a  ray  tracing  procedure 
is  used  to  calculate  ray  path  information  for  the  current  time  step;  this  information  is 
used  to  determine  the  expected  measurement  values  based  on  the  last  time  step. 
Second,  the  expected  measurement  values  are  subtracted  from  the  actual 
measurement  values.  This  information,  along  with  the  ray  path  geometry,  is  then 
preprocessed  so  that  it  can  be  used  by  the  neural  network.  Third,  the  neural  network  is 
trained  using  the  pre-  processed  data.  Fourth,  the  network  is  used  to  produce  an 
electron  density  estimate  for  the  current  time  step. 

Preliminary  experimental  results  indicate  that  ionospheric  electron  density  estimation 
using  this  method  is  both  computationally  feasible  in  real  time  and  capable  of  producing 
reasonable  reconstructions  given  the  proper  choice  of  neural  network.  In  particular,  a 
two  layer  neural  network  trained  used  the  Levenberg-Marquardt  method  was  found  to 
be  capable  of  locating  a  localized  time-varying  Gaussian  depletion  type  ionospheric 
disturbance  based  only  on  very  sparse  data. 

Second,  we  have  worked  on  the  development  of  a  new  fully  3D  reconstruction  algorithm 
for  CIT,  capable  of  combining  some  of  the  available  measurement  sources  in  order  to 
compensate  for  the  unavoidable  physical  limitations  of  the  acquisition  system. 

The  main  features  of  the  algorithm  are: 

a)  Fully  30  reconstruction  with  aititrary  voxel  size,  b)  Real-time  system  geometry 
computation  in  spherical  coordinates,  c)  Direct  (non-iterative)  reconstruction  based  on  a 
least-squares  inversion  method,  d)  Physical  constraints  (e.g.  non-negativity  of  the 
electron  density)  and  RADAR  measurements  (e.g.  WSBI  data)  introduced  as  boundary 
conditions  for  the  least-squares  algorithm,  e)  No  a  priori  information  on  the  electron 
density  distribution  required,  f)  Ready  for  acquisition  system  extension  to  occupation 
measurerrients.  When  satellite  occupation  TEC  measurements  will  become  available  an 
interface  is  ready  to  include  them  in  the  reconstruction  algorithm  system  matrix,  g) 
Choice  of  parameters  is  minimized.  Raw  estimates  of  the  measurement  errors  are  the 
only  parameters  necessary  to  run  the  whole  algorithm.  At  this  time  most  of  the 
software  features  are  operative.  The  package  includes  a  fast  3D  ray-tracing  algorithm 
capable  of  computing  a  system  matrix  with  arbitrary  voxel  sizes  in  real  time  and  storing 
it  in  a  compressed  sparse  form. 


The  bounded  least  squares  problem  engine  has  been  tested  and  optimized  for  the  given 
problem.  A  technique  to  improve  the  sparseness  of  the  intermediate  results  generated 
in  the  inversion  process  is  under  construction.  This  represents  at  the  moment  the  factor 
that  limits  the  size  of  the  system  matrix,  which  can  be  used  without  incurring  in 
unacceptable  memory  requirements. 

The  inclusion  of  non-TEC  measurement  data,  e.g.  WSBI  data,  is  under  development  in 
this  stage.  On  the  other  side  the  algorithm  gave  very  promising  results  simply  adding 
some  simple  physical  constraints  in  the  reconstruction  algorithm. 

Third,  we  developed  the  multi-input  volumetric  inversion  algorithm  (MIVIA).  After 
algorithm  development,  the  goals  of  the  research  were  to  test  its  reliability  and 
performance  for  various  types  of  CIT  system  geometry  and  various  ranges  of  regulatory 
parameters. 

Sensitivity  Analysis  of  MIVIA  involved  a  three  phase  sequence  of  tests  for  determining 
parametric  sensitivity,  convergence  sensitivity  and  sensitivity  to  a  priori  information. 

The  volumetric  inversion  algorithm  was  designed  to  select  a  priori  information  based  on 
guidance  from  TEC  measurements  obtained  during  the  CIT  satellite  overpass. 
Volurnetric  information  is  supplied  to  the  reconstruction  in  the  projection  and  image 
domains  of  each  iteration.  Investigation  of  the  relative  strengths  of  projection  and  image 
domain  corrections  yielded  optimal  values  of  1  for  projection  domain  shaping  and  0.01 
times  the  order  of  maximum  electron  density  for  image  domain  shaping.  Consequences 
of  excessive  image  domain  shaping  in  instances  of  low  sky  coverage  were  also 
demonstrated  including  the  case  of  over  correction  that  results  in  the  inversion  of  image 
features.  Convergence  was  achieved  faster  in  the  vertical  direction  than  along  the 
horizontal  because  of  the  dual  support  of  ray  path  profiles  and  image  neighborhood 
shapes  in  the  former.  Horizontal  convergence  supported  only  by  image  domain  shaping 
and  is  typically  achieved  between  40-50  iterations  as  opposed  to  10-15  iterations  for 
vertical  convergence.  With  respect  to  a  priori  input,  it  was  observed  that  reconstructions 
preserved  features  of  the  parent  image  when  aided  by  relatively  smooth  and  featureless 
a  priori  images.  Such  supportive  information  is  crucial  to  the  formation  of  the  vertical 
profile  and  is  used  based  on  the  need  for  such  Information  during  the  localized 
volumetric  shaping  processes. 

A  considerable  amount  of  work  was  performed  on  estimating  reconstruction  artifacts 
that  arise  because  of  non-ideal  system  geometry  and  incomplete  recording  of  an 
overpass  from  individual  stations.  Various  pathological  instances  of  imperfect  system 
geometry  were  developed  including  isolated  non-coplanar  receivers,  satellite  orbits 
longitudinally  displaced  from  receiver  chains  and  of  narrow  ray  path  beams  from  certain 
stations.  Among  the  important  observations  on  artifact  formation  was  that  displacement 
of  an  individual  station  from  an  otherwise  coplanar  receiver  chain  results  in  an  elevation 
of  features  in  the  reconstruction  in  over  the  latitudes  in  which  the  displaced  stations  ray 
fan  plays  a  dominant  role.  Artifacts  such  as  false  enhancements  and  troughs  were 
found  to  form  because  of  abrupt  termination  of  ray  fans  particularly  in  the  center  of  the 
image  plane.  It  was  also  observed  that  if  geographic  logistics  prohibit  the  placement  of 


a  receiver  chain  directly  below  the  satellite  path,  the  quality  of  reconstructions  is 
significantly  improved  by  positioning  at  least  one  additional  receiver  in  the  orbital  plane. 
Such  a  co-planar  receiver  provides  invaluable  support  to  the  set  of  largely  slant  TEC  ray 
paths  in  the  projection  set. 

A  scheme  for  the  quantitative  evaluation  of  the  features  of  an  ionospheric  image  were 
also  developed  with  the  aim  of  providing  an  objective  and  automatic  analysis  of  image 
content.  Histogram  redistribution  techniques  were  developed  to  separate  a  typical 
ionospheric  electron  distribution  into  enhancements  and  horizontally  stratified  bands.  In 
addition  extremely  localized  searches  were  devised  to  identify  travelling  ionospheric 
disturbances  which  are  usually  invisible  to  the  human  eye  because  of  the  dynamic 
range  of  the  entire  image.  Geometric  and  electron  density  based  properties  of  image 
regions  corresponding  to  features  such  as  peaks  and  horizontal  bands  were  derived 
and  used  for  the  comparison  of  images.  It  is  anticipated  that  this  quantitative  scheme 
will  permit  further  advanced  evaluation  of  CIT  inversion  techniques  by  comparing 
reconstructions  by  various  algorithms  using  the  same  TEC  data  set  or  for  comparing 
CIT  reconstructions  to  gold  standards  such  as  images  obtained  using  the  incoherent 
scatter  radar  technology. 

Finally,  we  worked  on  developing  and  implementing  CVT  (CREDO-VSIRT  Technique), 
studies  on  the  incorporation  of  WSBI  data  into  VSIRT,  and  reconstructions  for  simulated 
and  real  data. 

The  studies  on  the  incorporation  of  WSBI  data  into  VSIRT  entailed  understanding  the 
weaknesses  of  both  CREDO  and  VSIRT,  and  finding  ways  to  use  the  strengths  of  each 
for  the  benefit  of  better  reconstructions.  Some  of  the  techniques  include  the  initial 
guess  method,  the  control  point  method  and  the  method  of  doing  the  WSBI  inversion 
within  VSIRT. 

In  the  area  of  reconstructions,  the  most  important  result  was  the  simulated  test  of 
detecting  an  E-layer  using  VSIRT.  This  test  proved  that  VSIRT  was  able  to  pick  up  an 
E-layer,  but  it  was  only  able  to  do  it  in  a  very  controlled  environment,  i.e.  ideal  station 
alignment  and  an  E-layer  that  was  quite  high  in  altitude. 

Technical  details  are  provided  in  the  enclosed  papers. 
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Short  title:  ESTIMATING  IONOSPHERIC  CHANGES 


Abstract.  We  present  a  novel  artificial  neural  network  based  computerized 
ionospheric  tomography  (CIT)  technique,  capable  of  imaging  through  time  in  three 
dimensional  space.  Total  electron  content  (TEC)  data  collected  from  a  satellite  passing 
over  the  region  of  interest  is  used  to  train  a  neural  network.  The  trained  network 
creates  estimates  of  the  difference  in  ionospheric  electron  density  between  time  steps 
based  on  the  TEC  data.  Application  of  the  difference  estimate  at  each  time  step  to  the 
previous  electron  density  image  results  in  a  time  varying  ionospheric  electron  density 
estimate.  Experimental  results  on  synthetic  data  are  presented  that  demonstrate  that 
the  algorithm  is  capable  of  detecting  short-term  localized  disturbances  in  the  ionosphere. 


1.  Introduction 


Although  accurate  information  about  ionospheric  electron  densities  is  crucial  for 
correcting  radio  signal  aberrations,  the  ionosphere  is  far  from  static,  varying  even  within 
the  time  scale  of  a  single  satellite  pass  over  the  region  of  interest.  Natural  effects  such 
as  diurnal  variation  along  the  day/night  boundary,  solar  radio  emissions,  solar  flares, 
volcanic  eruptions,  earthquakes,  and  geomagnetic  storms  can  all  affect  the  dynamics  of 
the  ionosphere  within  minutes.  Rapid  perturbation  of  the  ionosphere  can  also  result 
from  man  made  effects  such  as  chemical  and  nuclear  explosions  and  high-powered  radio 
waves. 

Current  measurement  techniques  are  unable  to  provide  the  necessary  information 
about  ionospheric  electron  densities.  As  a  result,  it  is  necessary  to  use  computational 
methods  to  transform  measurement  data  into  images  of  ionospheric  electron  densities. 
The  use  of  such  computational  techniques  is  referred  to  as  computerized  ionospheric 
tomography  (CIT). 

The  majority  of  current  CIT  algorithms  are  only  capable  of  reconstructing  static 
images  of  the  ionosphere.  However,  when  static  CIT  algorithms  are  applied  to  a 
dynamic  ionosphere,  the  clash  between  the  assumption  of  time  invariance  made  by 
these  algorithms  and  the  time  dependency  of  the  actual  electron  density  results  in 
estimation  errors.  These  errors  include  clearly  visible  artifacts,  such  as  streaking, 
negative  ionospheric  density  values,  and  geometric  distortion  (i.e.,  features  displaced 
from  their  actual  position)  [Sutton  and  Na,  1998]. 

In  order  to  avoid  such  errors,  it  is  necessary  to  use  time  varying  CIT  (TVCIT) 
algorithms  capable  of  estimating  ionospheric  electron  densities  as  they  evolve  in  time. 
Methods  have  been  presented  that  attempt  to  estimate  time  varying  ionospheric  images 
[Sutton  and  Na,  1998;  Howe  et  al.,  1998].  However,  these  methods  are  critically 
dependent  on  a  priori  models  of  ionospheric  structure.  As  a  result,  ionospheric  motion 
and  structure  that  does  not  fit  the  a  priori  models  may  be  ignored  or  misinterpreted. 


Unfortunately,  it  is  exactly  this  kind  of  unpredictable  or  unknown  behavior  that  it  is 
most  desirable  to  detect. 

This  problem  stems  ultimately  from  the  underdetermined  nature  of  the  system. 
The  above  methods  attempt  to  address  the  lack  of  data  by  making  use  of  empirically 
derived  basis  functions.  However,  in  doing  so  they  compromise  their  ability  to  detect 
novel  or  unusual  ionospheric  structure  and  activity.  In  this  paper  we  present  an 
alternative  artificial  neural  network  based  technique  capable  of  imaging  through  time  in 
three  dimensional  space.  The  neural  network  based  method  attempts  to  preserve  the 
ability  to  detect  interesting  ionospheric  phenomenon,  while  avoiding  the  problem  of  data 
scarcity.  Because  the  data  is  so  sparse,  it  is  extremely  difficult  to  accurately  compute 
the  electron  density  at  any  point  in  the  image.  Instead,  we  seek  to  generate  estimates 
of  the  ion  density  based  on  what  little  information  is  available.  In  principle,  artificial 
neural  networks  are  well  suited  to  this  kind  of  estimation  because  they  have  the  ability 
to  generalize  results  based  on  a  small  subset  of  the  data. 

2.  Computerized  Ionospheric  Tomography 

A  minimum  data  collection  system  used  for  CIT  consists  of  a  radio  beacon  satellite 
orbiting  over  the  region  of  interest  and  a  series  of  ground  station  receivers  aligned 
(roughly)  along  the  satellite  path.  The  system  collects  measurements  of  total  electron 
content  (TEC)  by  means  of  the  Doppler  radar  technique.  Two  signals  are  sent  out  at 
frequencies  high  enough  to  linearly  penetrate  the  ionosphere;  the  phase  difference  at 
each  ground  station  receiver  provides  a  measure  of  the  TEC  along  the  ray  path  from 
the  satellite  to  that  station.  This  TEC  measurement  is  approximately  the  line  integral 
of  electron  densities  along  the  ray  path. 

Figure  1 

The  first  CIT  reconstruction  technique  was  introduced  by  Austen  et  al.,  1986.  Since 
that  time,  a  number  of  CIT  techniques  have  been  devised  and  experiments  in  various 
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locations  around  the  world  have  been  conducted. 


3.  Neural  Network  Approach  to  TVCIT 

The  artificial  neural  network  is  not  used  to  estimate  the  electron  density  directly 
but  rather  the  change  in  the  electron  densities  at  each  time  step.  If  the  time  step  used 
is  sufficiently  small,  then  the  difference  in  ion  density  between  time  steps  should  be 
minor.  Thus,  the  problem  can  be  broken  down  into  a  series  of  relatively  small  estimation 
problems. 

At  each  point  in  time,  the  algorithm  proceeds  through  four  main  steps.  First,  a 
ray  tracing  procedure  is  used  to  calculate  ray  path  information  for  the  current  time 
step;  this  information  is  used  to  determine  the  expected  TEC  values  based  on  the  last 
time  step.  Second,  expected  TEC  values  are  subtracted  from  the  actual  measured  TEC 
values  for  the  current  time  step;  this  information,  along  with  the  ray  path  geometry,  is 
preprocessed  so  that  it  can  be  used  by  the  neural  network.  Third,  the  neural  network  is 
trained  using  the  preprocessed  data.  Fourth,  the  network  is  used  to  create  the  electron 
density  estimate  for  the  current  time  step. 

The  algorithm  assumes  that  a  good  initial  estimate  is  available  at  the  start  of  the 
procedure.  Although  this  is  a  significant  assumption,  there  are  a  few  practical  methods 
by  which  this  initial  estimate  could  be  obtained;  these  will  be  discussed  later. 

Aside  from  an  initial  estimate,  the  algorithm  also  requires  system  geometry 
information  and  TEC  measurement  data.  The  geometry  is  determined  by  the  satellite 
position  and  ground  stations  positions.  The  TEC  measurement  data  is  collected  by  the 
ground  stations  at  select  time  intervals. 

4.  Algorithm  Details 

Given  the  system  geometry  and  TEC  data,  the  first  step  in  the  reconstruction 
process  is  the  calculation  of  the  ray  path  information.  In  general  in  tomography 


problems  one  seeks  to  calculate  a  system  matrix  (either  explicitly  or  implicitly)  which 
contains  information  about  the  ray  paths  through  the  region  being  imaged.  In  the 
ionospheric  tomography  case,  at  each  instant  in  time  there  is  a  spread  of  ray  paths  from 
the  satellite  to  each  of  the  ground  stations.  The  data  used  is  collected  from  within 
a  given  time  window.  Once  that  time  window  has  been  selected,  it  is  necessary  to 
calculate  a  system  matrix  describing  the  ray  paths  from  each  of  the  satellite  positions  to 
the  ground  stations.  Calculating  the  system  matrix  in  general  involves  some  form  of  ray 
tracing.  In  this  case,  the  ray  tracing  procedure  is  somewhat  involved  due  to  the  nature 
of  the  spherical  coordinate  system  used  by  the  simulation. 

The  ray  path  information  is  then  used  to  calculate  expected  TEC  values  given  the 
electron  density  distribution  of  the  previous  time  step.  These  expected  values  are  then 
subtracted  from  the  measured  TEC  values  to  produce  TEC  difference  values. 

The  next  step  of  the  algorithm  consists  of  preprocessing  of  the  ray  path  and 
measurement  data.  Preprocessing  of  the  data  is  necessary  because  the  artificial 
neural  network  is  expected  to  approximate  the  change  in  ion  density  as  a  function  of 
spatial  position.  As  such,  it  requires  as  its  input  pixel  coordinates  and  outputs  the 
corresponding  electron  density  change  value.  This  requires  that  the  data  (i.e.,  the  total 
electron  content  measurement  differences  and  system  matrix)  be  used  to  generate  pixel 
values  on  which  to  train  the  artificial  neural  network.  This  is  done  by  calculating  the 
the  minimum  norm  solution  to  the  underdetermined  equation 

Ax  =  b, 

where  A  is  the  sparse  system  matrix,  b  is  a  vector  containing  the  TEC  difference 
measurements,  and  x  is  a  vectorized  representation  of  the  image.  The  resulting  solution 
vector  X  contains  the  training  output  values,  each  of  which  corresponds  to  a  non-zero 
element  of  the  system  matrix.  The  coordinates  of  the  non-zero  elements  are  the  training 
input  values.  Note  that  these  are  only  the  pixel  values  along  the  ray  paths;  the  procedure 


is  not  capable  of  calculating  of  calculating  pixel  values  elsewhere  in  the  image. 

Once  the  experimental  data  has  been  preprocessed  into  the  correct  form,  it  can  be 
used  to  train  the  artificial  neural  network.  The  network  has  three  inputs,  one  for  each 
of  the  coordinate  values,  and  a  single  output  unit  for  the  electron  density  change.  The 
output  unit  has  a  linear  response  so  that  it  is  capable  of  representing  any  value.  The 
number  and  nature  of  the  hidden  units  varies  by  the  approach  used. 

The  network  is  then  trained  using  the  Levenberg-Marquardt  method.  An  alternative 
to  the  gradient  descent  based  backpropagation  method,  the  Levenberg-  Marquardt 
algorithm  is  theoretically  faster  and  potentially  more  accurate,  although  it  requires 
more  storage  space  for  the  computation. 

The  Levenberg-Marquardt  algorithm  can  be  used  as  an  approximation  to  the 
Newton  method.  The  weight  update  rule  is 

AW  =  (J^J -f /iI)"^J^e, 

where  W  is  a  matrix  of  weights  Wij  connecting  two  layers,  J  is  the  Jacobian  matrix  of 
derivatives  of  each  error  to  each  weight,  e  is  an  error  vector,  and  controls  the  learning 
rate. 

If  fj.  is  very  large,  the  update  approximates  gradient  descent;  if  it  is  small  the 
update  approximates  the  pure  Newton  method.  The  Newton  method  approximation  is 
faster,  but  tends  to  be  less  accurate  when  near  the  minima.  To  deal  with  this  problem, 
the  parameter  n  is  adjusted  in  an  adaptive  fashion.  If  the  error  is  decreasing,  n  is  made 
bigger;  if  the  error  is  increasing,  n  is  made  smaller. 

When  training  is  complete,  a  complete  estimated  ion  density  change  image  can  then 
be  created.  The  coordinates  of  each  pixel  are  given  as  input  to  the  net  and  the  output  is 
assigned  as  the  pixel  value,  thus  forming  a  complete  image.  The  electron  density  change 
image  can  then  be  added  to  the  previous  electron  density  image  to  produce  the  current 
electron  density  image. 


5.  Experimental  Results 


Experimental  results  for  the  algorithm  were  generated  using  simulated  volume  data 
along  with  ground  station  positions  and  satellite  paths  from  an  ongoing  campaign  in 
the  Caribbean  region. 

The  simulation  included  an  ionospheric  disturbance,  a  Gaussian  depletion  of  the 
electron  density.  Initially,  there  was  no  activity.  Then  a  region  of  increased  electron 
density  was  introduced,  followed  soon  after  by  a  neighboring  region  of  decreased  electron 
density. 

Figures  X.X  and  X.X  show  a  longitudinal  slice  of  the  electron  density  image 
produced  by  the  simulation  before  the  disturbance. 

Figure  2  figure  3 

ionospheric  density  estimation  problem  appears  to  yield  good  localization  of  an 
ionospheric  disturbance,  as  long  as  some  of  the  ray  paths  pass  close  enough  to  the 
disturbance  for  it  to  affect  the  measurement  data. 

Figure  X.X  shows  the  output  of  the  Levenberg-Marquardt  trained  network  with 
ten  neurons  given  four  ray  paths,  one  of  which  passes  close  enough  to  the  Gaussian 
depletion  to  detect  it  at  that  time  step.  Figure  4.15  should  be  compared  to  figure 
X.X,  which  displays  the  ground  truth  electron  density  change  for  the  same  time  and 
longitudinal  slice.  Although  the  electron  density  change  is  more  distributed  than  in  the 
ground  truth,  the  peak  of  the  density  change  is  located  in  the  correct  spot. 

Figure  4  Figure  5 

Figures  X.X-X  show  the  output  of  that  same  network  at  the  1st,  7th,  and  16th 
longitudinal  slices,  along  with  their  associated  ground  truths.  Note  that  the  electron 
density  is  located  correctly  in  longitude  as  well  as  in  latitude  and  altitude. 

Figure  6  Figure  7  Figure  8  Figure  9  Figure  10  Figure  11 

The  method  suffers  from  the  fact  that  the  number  of  hidden  units  must  be  chosen 
carefully.  If  too  few  hidden  units  are  used,  the  resolution  is  insufficient;  if  too  many 
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hidden  units  are  used  the  network  overfits  the  data  and  estimate  becomes  distorted  and 
irregular. 

As  well,  good  convergence  via  the  Levenberg-Marquardt  method  is  sensitive  the 
random  initialization  of  the  neural  network  weights.  Therefore,  for  the  ionospheric 
reconstruction  system  to  produce  reliable  results,  it  is  necessary  to  repeatedly  train  the 
network  at  each  time  step.  The  result  at  a  given  time  step  with  the  least  error  is  then 
used.  Although  this  still  does  not  guarantee  a  good  approximation  (it  might  be  that  a 
poor  choice  of  random  starting  weights  is  made  each  each  try),  it  would  greatly  reduce 
the  probability  of  a  poor  approximation.  The  number  of  tries  can  be  adjusted  according 
the  required  confidence  in  the  results. 

However,  such  an  approach  increases  the  amount  of  computation.  Fortunately, 
training  the  network  using  the  Levenberg-Marquardt  method  takes  considerably  less 
time  than  generating  the  complete  ion  density  estimate  from  the  trained  network.  On 
the  order  of  ten  training  attempts  can  be  made  without  adversely  affecting  the  speed  of 
the  algorithm  to  a  significant  degree. 

6.  Conclusions  and  Future  Work 

We  have  presented  a  novel  neural  network  based  CIT  method  capable  of  imaging 
ionospheric  electron  densities  across  time  and  three  dimensional  space.  In  experiments 
using  simulated  electron  density  data  and  realistic  CIT  system  geometries,  the  algorithm 
has  proven  capable  of  estimating  the  location  and  magnitude  of  regional  transient 
disturbances  in  ionospheric  electron  density.  The  method  has  the  advantage  of  not 
relying  on  a  priori  models  of  the  ionosphere,  thereby  giving  it  the  capability  to  detect 
unusual  or  unexpected  ionospheric  phenomenon.  Furthermore,  the  computation  can 
be  performed  in  real  time,  making  it  amenable  to  field  application.  The  method  is 
extenable  as  well  in  that  it  can  be  easily  modified  to  use  additional  data  sources,  such  as 
ionosonde  readings  or  total  electron  content  measurements  garnered  from  GPS  satellites. 


Still,  some  improvements  suggest  themselves  with  regard  to  this  TVCIT  method. 
The  sensitivity  of  the  method  to  the  random  initialization  can  be  mitigated  by  developing 
a  system  which  attempts  to  pick  the  best  result  from  an  array  of  estimates  generated 
by  different  initializations.  Or  alternatively,  if  the  time  resolution  is  sufficiently  fine, 
then  the  electron  density  change  can  be  considered  continuous  in  time,  and  information 
about  the  motion  in  previous  times  could  be  used  to  initialize  the  network  weights  for 
the  current  time  step. 

As  well,  possible  overfitting  of  the  neural  network  to  the  limited  data  could  be 
reduced.  This  could  be  done  by  sampling  the  data  according  to  ray  path  coverage 
density,  or  by  taking  advantage  of  the  extendability  of  the  method  and  adding  additional 
data  sources. 

Effective  initialization  of  the  algorithm  remains  the  most  pressing  problem.  In  the 
experiments  discussed  previously,  a  good  initial  estimate  was  assumed.  However,  a  poor 
estimate  can  ruin  the  results  of  the  algorithm.  Incoherent  scattering  radar  measurements 
could  provide  a  good  initial  estimate;  however  in  the  absence  of  such  generally  available 
measurements  it  is  necessary  to  develop  a  more  general  initialization  method. 

Three  such  methods  suggest  themselves.  The  first  would  be  to  start  with  a  set  of 
a  priori  images,  chosen  appropriately  for  the  geographical  placement  of  the  region  of 
interest.  The  method  could  then  be  initialized  with  the  a  priori  image  that  most  closely 
matches  the  measurement  data  from  the  first  time  step. 

The  second  approach  would  be  to  collect  measurement  data  for  some  amount  of 
time  as  the  satellite  begins  its  pass  over  the  region  of  interest.  A  standard  non-time 
varying  CIT  technique  could  then  be  used  to  create  a  reconstruction  using  the  data 
collected  from  the  beginning  of  the  pass;  the  resulting  image  could  then  be  used  to 
initialize  the  algorithm. 

The  third  method  would  only  be  applicable  in  a  scenario  where  repeated  experiments 
are  conducted.  During  the  first  satellite  pass,  only  a  static  reconstruction  would  be 


performed.  This  static  reconstruction  would  then  be  used  to  initialize  the  time  varying 
algorithm  for  the  next  satellite  pass. 
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Figure  1.  A  simple  CIT  system. 

Figure  2.  A  2-dimensional  view  of  the  10th  longitudinal  slice  of  the  simulated  electron 
density. 

Figure  3.  A  3-dimensional  view  of  the  10th  longitudinal  slice  of  the  simulated  electron 
density. 

Figure  4.  10th  longitudinal  slice  of  the  electron  density  change  estimate. 

Figure  5.  10th  longitudinal  slice  of  the  simulated  electron  density  change. 

Figure  6.  First  longitudinal  slice  of  the  electron  density  change  estimate. 

Figure  7.  First  longitudinal  slice  of  the  simulated  electron  density  change. 

Figure  8.  7th  longitudinal  slice  of  the  electron  density  change  estimate. 

Figure  9.  7th  longitudinal  slice  of  the  simulated  electron  density  change. 

Figure  10.  16th  longitudinal  slice  of  the  electron  density  change  estimate. 

Figure  11.  16th  longitudinal  slice  of  the  simulated  electron  density  change. 
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Abstract 

This  paper  presents  results  of  investigations  on  the  relationship  between  optimal  image 
resolution  and  associated  sky-coverage  for  ionospheric  tomography.  Several  measures  of  the 
variation  of  image  information  content  as  a  function  of  resolution  have  been  presented  for 
understanding  the  nature  of  tradeoffs  between  information  content  and  computational  resources. 
A  new  quantitative  description  of  imaging  geometry  enables  sky-coverage  analysis  for  this 
limited  information  tomographic  system.  The  paper  also  presents  examples  of  commonly 
encountered  non-ideal  system  geometries,  diagnosis  of  associated  image  artifacts,  and 
improvements  achieved  by  minor  logistic  modifications  to  the  imaging  system. 
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1.  Introduction 

Computerized  ionospheric  tomography  (CIT)  [Austen  et  al,  1988]  provides  an  excellent 
opportunity  for  observing  the  relationship  between  the  information  content  registered  by  a 
limited  angle,  limited  data  tomographic  imaging  system  and  the  subsequent  information 
distribution  achieved  in  the  tomographic  inversion  process.  The  availability  of  records  from 
numerous  experimental  campaigns  conducted  over  the  past  decade  and  the  emergence  of 
applications  such  as  over-the-horizon  radar  provide  additional  incentive  for  the  study  of 
information  flow  in  restricted  tomographic  systems.  Knowledge  of  information  availability  and 
distribution  in  limited  information  systems  is  expected  to  contribute  significantly  to  diverse 
pursuits  ranging  from  optimization  of  computational  resources  to  the  diagnosis  of  misleading 
features  of  artifacts  in  the  reconstructions.  This  paper  is  concerned  with  two  specific  sets  of  CIT 
related  analyses  governed  by  tomographic  information  content,  namely,  image  resolution 
analysis  and  sky-coverage  based  image  artifact  analysis. 

Ionospheric  tomography  involves  registration  of  the  phase  difference  between  two  radio 
signals  from  a  low-orbit  satellite  by  one  or  more  chain  of  ground  based  receivers.  The  phase 
difference  being  representative  of  the  total  electron  content  (TEC)  along  the  corresponding  signal 
ray-path  through  the  atmosphere,  the  recorded  data  can  be  processed  using  tomographic 
principles  to  reconstruct  a  two  or  three-dimensional  image  of  the  ionospheric  region  lying 
between  the  satellite  and  the  ground  stations.  Due  to  geometric  restrictions  such  as  the  curvature 
of  the  earth  and  continuity  of  the  atmospheric  shell  around  it  and  logistic  restrictions  such  as  the 
limited  availability  and  irregular,  non-coplanar  positioning  of  receivers,  CIT  is  effectively  a 
limited  angle,  limited  sample  tomographic  system.  In  addition  to  the  above  mentioned  static 
restrictions,  variations  in  receiver  operation  often  result  in  dynamic  limitations  such  as  loss  of 
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absolute  phase  related  information  [Leitinger,  1994]  and  the  incomplete  recording  of  satellite 
overpasses  from  individual  stations.  Therefore  the  comprehensive  study  of  information 
registration  and  distribution  would  involve  quantitative  definitions  of  ray-path  distribution  or 
sky-coverage  in  the  ionospheric  region  of  interest.  The  optimal  resolution  of  reconstructed 
images  in  limited  angle  and  limited  sample  imaging  systems  is  governed  by  sky-coverage.  In 
addition,  the  distribution  of  sky  coverage  itself  can  be  related  to  the  formation  of  image  artifacts 
because  it  is  effectively  a  measure  of  limited  visibility  of  the  object  distribution. 

The  remainder  of  this  paper  is  organized  into  four  sections.  A  more  detailed  description 
of  sky-coverage  and  the  image  coordinate  system  are  presented  in  section  two.  Sections  three 
and  four  discuss  the  important  results  and  observations  made  in  the  process  of  image  resolution 
and  image  artifact  analysis  based  on  sky-coverage.  The  concluding  section  consists  of 
discussions  on  the  significance  of  information  content  evaluation  and  artifact  diagnosis  for 
limited  angle  tomographic  systems  and  on  directions  for  future  research  in  this  area. 

2.  The  Coordinate  System,  Coverage  and  Imaging  Geometry 

Traditionally  CIT  systems  have  used  satellites  with  polar  orbits  and  therefore  ground 
based  receiver  chains  have  usually  been  oriented  in  a  roughly  north-south  direction.  An 
exhaustive  study  of  experimental  imaging  geometries  used  over  the  past  decade  shows  major 
digressions  from  the  assumptions  of  chain  linearity  and  coplanarity  of  ground  stations  and  the 
orbital  plane  [Biswas,  1999].  Apart  from  irregular  latitudinal  spacing  of  receivers  and 
longitudinal  displacement  of  individual  stations  from  the  orbit  plane,  there  exists  a  wide  variation 
of  ground  chain  distributions  and  alignments,  the  usage  of  multiple  chains  and  the  inclusion  of 
multiple  satellite  orbits  in  the  imaging  process.  Further  more,  it  is  necessary  to  make  provisions 
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for  the  usage  of  non-polar  satellites  in  future  systems.  The  practice  of  converting  non-coplanar  or 
slant  ray-paths  to  vertical  coplanar  ray  paths  together  with  adjustments  of  TEC  values  involves 
geometric  and  electron  density  related  assumptions  and  extrapolations  [Mitchell  et  al,  1997]. 
Recent  multiple  chain  CIT  campaigns  have  used  numerical  interpolations  between  a  set  of 
verticalized  two-dimensional  image  planes  to  obtain  a  three-dimensional  reconstruction. 
However,  as  demonstrated  by  a  volumetric  inversion  technique  for  CIT,  it  is  most  convenient  to 
use  the  true  geocentric  three-dimensional  geometry  for  CIT  imaging  since  it  avoids  data 
preprocessing  and  permits  direct  use  of  slant  TEC  values.  The  three  coordinates  of  the  geocentric 
coordinate  system  are  latitude  (X-axis),  altitude  (Y-axis)  and  longitude  (Z-axis).  Whereas  the 
resolution  can  be  conveniently  varied  to  desired  resolutions  with  uniformity  for  the  first  two 
coordinates,  the  latitudinal  dependence  of  longitudinal  spacing  makes  it  marginally  more 
computationally  intensive  to  ensure  uniform  spacing  in  the  Z  direction.  Fortunately,  fairly  low 
resolution  of  between  3-7  planes  is  sufficient  for  most  CIT  imaging  systems. 

For  each  satellite  overpass,  there  exists  a  unique  sky-coverage  map.  Coverage  is  defined 
as  an  image,  occupying  the  same  geographic  region  as  the  tomographic  image,  whose  voxel 
value  is  a  count  of  the  number  of  ray-paths  passing  through  it.  Thus  the  coverage  map  shows  the 
contribution  of  information  in  the  most  primary  sense  from  various  portions  of  the  ionosphere  to 
the  TEC  data  set.  Figure  1  shows  a  diagrammatic  comparison  between  the  imaging  systems  and 
coverage  maps  of  typical  medical  and  ionospheric  systems.  The  limitations  of  angular  visibility, 
number  and  regularity  of  samples  are  clearly  evident  in  the  second  imaging  system.  The 
rectangular  image  of  Figure  Id  shows  curved  ray-paths  because  of  the  curved  to  flat  earth 
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transformation. 
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The  system  devised  for  the  quantitative  description  of  CIT  imaging  geometry  is 
represented  in  Figure  2.  The  suite  of  coverage  parameters,  enlisted  and  explained  in  Table  1,  can 
be  classified  into  satellite  orbit  and  ground  station  parameters.  These  parameters  permit 
characterization  and  comparison  of  CIT  imaging  systems.  They  have  been  used  in  designing 
several  simulated  geometries  for  coverage  analysis  in  section  four. 

3.  Resolution  Analysis  for  Ionospheric  Tomography 

The  optimal  degree  of  three-dimensionality  for  volumetric  reconstruction  depends  on  the 
tradeoff  between  the  feature  resolution  and  information  density,  that  is,  coverage.  TEC  records 
from  four  CIT  campaigns  were  analyzed  for  coverage  estimates  while  varying  voxel  resolution 
along  various  image  coordinates  for  a  fixed  geographic  volume.  All  four  campaigns  were 
conducted  in  the  equatorial  and  mid-latitudes,  hence,  there  was  significant  gain  in  earth-surface 
resolution  with  small  increments  in  the  number  of  longitudinal  planes.  Figure  3  shows  the  system 
geometries  of  the  four  campaigns  which  represent  the  geometric  diversity  encountered  in  CIT. 

As  the  resolution  was  varied  along  the  three  geocentric  coordinates  with  geographic 
resolutions  shown  in  Table  2,  there  was  an  initial  decrease  in  coverage  with  increase  in 
resolution  followed  by  convergence  as  demonstrated  in  Figure  4.  The  point  of  coverage 
stabilization  was  regarded  as  the  optimal  resolution  along  the  corresponding  coordinate  axis  for 
the  given  CIT  geometry. 

The  results  of  resolution  analysis  have  been  presented  in  Table  3  with  Figure  5 
highlighting  some  key  observations.  Comparing  coverage  of  various  image  resolutions  along  the 
latitudinal  and  altitudinal  directions  for  the  simple,  single  chain  imaging  systems  of  the  Mid¬ 
American  err  Experiment  of  1993  (MACE-93)  and  the  Australian  campaign  of  1995,  it  is 
observed  that  optimal  resolution  along  the  latitudinal  direction  is  governed  largely  by  the 
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displacement  of  the  satellite  orbit  from  the  ground  station  and  in  the  altitudinal  direction  by  the 
density  and  linearity  of  the  receiver  chain.  MACE-93  presented  a  fairly  dense  and  linear  receiver 
chain  but  one  that  was  longitudinally  displaced  from  the  satellite  orbit  whereas  the  reverse  was 
the  case  with  the  Australian  system.  Consequently,  there  is  gainful  return  for  expending 
computational  resources  towards  improving  the  resolution  along  the  X  axis  for  the  former  and 
along  the  Y  axis  for  the  latter.  Data  from  Table  3  indicates  that  given  typical  TEC  ray-path 
density  of  a  CIT  system,  optimal  voxel  sizes  lie  between  1 1 1  and  150  km  along  the  X  axis  on  the 
earth  surface  and  around  20  km  along  the  Y  axis  or  altitudinal  direction. 

Using  a  latitudinal  resolution  of  approximately  120  km  and  altitudinal  resolution  of  20 
km.  Figure  6  shows  the  variation  of  percentage  coverage  along  each  of  the  longitudinal  or  Z 
planes  for  the  four  imaging  systems.  The  single  most  important  observation  is  that  even  for 
conventional  single  chain,  single  orbit  geometries,  such  as  MACE-93  and  Austraiia-95,  there  is 
significant  coverage  in  non-central  longitudinal  planes.  This  justifies  the  use  of  multiple 
longitudinal  planes  and  hence  true  three-dimensional  coordinates  for  CIT  inversion  particularly 
for  multiple  chain  geometries  such  as  for  the  ionospheric  correction  test  methods  II  system  of 
1996  (ICMT2-96)  and  of  multi -chain,  multi-orbit  geometries  such  as  used  in  the  Caribbean 
campaign  of  1997.  Another  noteworthy  aspect  of  Z-plane  information  is  that  coverage  is  highest 
under  orbital  planes  and  the  change  in  sub-orbital  coverage  is  affected  to  a  greater  extent  by 
orbit-chain  displacement  as  evidenced  by  smaller  changes  in  coverage  for  the  Australian 
campaign  than  for  MACE-93  for  Z  numbers  of  3  and  1 1 . 

Although  the  definition  of  coverage  used  in  this  study  did  not  account  for  the  angular 
visibility  of  image  voxels  and  only  the  percentage  coverage  was  used  for  resolution  analysis,  it  is 
an  estimate  of  the  absolute  maximum  information  content  of  the  measured  TEC  data.  This  study 
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revealed  the  significance  of  system  characteristics  such  as  chain  linearity  and  orbital 
displacement  in  determining  available  information  content  for  processing  and  therefore,  in 
serving  as  a  measure  of  optimal  image  resolution. 

4.  Coverage  Analysis  for  Ionospheric  Tomography 

A  study  of  several  recent  CIT  campaigns  {Biswas,  1999]  show  that  the  most  common 
limitations  of  CU  systems  are  limited  availability  and  irregular  spacing  of  receivers,  longitudinal 
displacement  of  individual  receivers  or  the  entire  station  chain  from  the  orbital  plane  and  the 
incomplete  recording  of  the  overpass  from  individual  stations  resulting  in  abrupt  edges  in  the 
coverage  map.  In  a  typical  CIT  system  intended  to  cover  40-50°  of  latitude,  approximately  1200 
project  samples  are  recorded  by  4-7  receivers  within  ±45°  of  the  vertical  axis  of  the  image.  This 
results  in  a  dense  overlap  of  ray-paths,  hence  high  coverage,  in  the  center  of  the  image  which 
gradually  declines  to  null  at  either  lateral  margin  of  the  image  providing  a  U-shaped  cone  of 
visibility  as  shown  in  Figure  1.  In  CIT,  the  typical  receiver  spacing  of  100-300  km  on  the  earth 
surface  is  significant  when  compared  with  source-receiver  separation  of  only  1000-1200  km. 
According  to  the  fundamentals  of  tomography  [Kak  and  Slaney,  1988],  such  severe 
undersampling  can  result  in  significant  reduction  in  high  frequency  energy  of  the  projection  set 
and  by  extension,  of  the  reconstructed  image.  Images  of  ionospheric  electron  content  are,  by 
nature,  smooth  distributions  of  plasma  and  of  low  spatial  frequency  content.  Hence,  visual 
separation  of  genuine  image  features  from  artifacts  due  to  non-ideal  sampling  for  CIT  poses  an 
additional  challenge.  Sky-coverage  maps  provide  useful  information  about  the  relative  visibility 
of  various  regions  of  the  ionosphere  to  the  imaging  system  during  a  satellite  overpass  and  can  be 
directly  related  to  the  formation  of  image  artifacts. 
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The  principle  aim  of  experiments  on  sky-coverage  for  ionospheric  tomography  was  the 
diagnosis  of  common  image  artifacts  resulting  from  non-ideal  aspects  of  the  imaging  geometry. 
As  discussed  in  the  first  two  sections,  the  study  was  based  on  devising  various  pathological 
situations  which  further  restrict  the  information  content  of  this  inherently  limited  angle,  limited 
data  imaging  system  of  CIT. 

Simple  geometric  images  approximating  prominent  ionospheric  formations  such  as  peaks 
embedded  in  horizontally  stratified  layers  with  or  without  horizontal  gradients  were  developed 
for  easy  diagnosis  and  appraisal  of  reconstruction  artifacts  resulting  from  the  above  mentioned 
limitations.  The  multi-source  volumetric  inversion  algorithm  (MIVIA)  [Biswas  and  Na,  1998a, 
b]  used  for  the  reconstruction  process  used  the  test  image  itself  to  draw  a  priori  information 
based  on  the  data  obtained  in  the  simulated  CIT  overpass  for  various  imaging  geometries. 
Volumetric  reconstruction  uses  shape  information  from  electron  distributions  of  the  a  priori 
images  to  guide  CIT  reconstruction. 

Figures  7a  through  7e  represent  the  simulated  CIT  imaging  systems  and  coverage  maps 
associated  with  the  data  sets  used  in  reconstructions  of  Figures  8b  through  8f.  Despite  the  close 
match  of  a  priori  and  simulated  TEC  data,  the  non-ideal  geometries  of  all  the  systems  except  that 
of  Figure  7a,  result  in  misleading  features  or  artifacts  in  the  reconstructions.  Figure  8a  shows  the 
original  test  image  and  Figure  8b,  the  reconstruction  obtained  with  a  near-ideal  five-station 
geometry. 

The  reconstructed  image  of  Figure  8c  corresponds  to  the  three-station  geometry  of  Figure 
7b  where  the  second  station  is  displaced  from  the  satellite  orbit.  The  coverage  from  the  two  sub¬ 
orbital  stations  is  incomplete  and  there  is  a  region  above  which  the  only  coverage  is  from  the 
second  station.  The  corresponding  region  of  the  reconstruction  shows  an  elevation  in  the  altitude 


Figure  7 


Figure  8 
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of  horizontal  band.  This  altitudinal  shift  can  be  explained  by  two  geometric  the  two  contributing 
factors  as  shown  in  Figure  9.  The  first  of  these  factors  is  the  increase  in  the  conditional  number 
or  the  ratio  of  the  largest  to  the  smallest  eigenvalue  of  the  sub-system  matrix  corresponding  to 
ray-paths  from  the  displaced  ground  stations.  As  shown  in  Figure  9a,  the  farther  a  receiver  is 
from  the  satellite  orbit,  the  narrower  the  angle  subtended  by  a  given  segment  of  the  orbit  on  the 
ground  and  the  smaller  the  range  of  angles  over  which  projection  sets  are  available  at  that  site  in 
the  tomographic  sense.  A  smaller  number  of  angles  given  the  same  number  of  TEC  ray  paths 
accounts  for  increase  in  the  overall  under-determinedness  of  the  inversion  problem  and  translates 
into  lack  of  tomographic  information  from  the  under-represented  region  of  the  atmosphere.  The 
second  reason  for  the  altitude  shift  can  be  attributed  to  the  discretization  of  the  ray-path  in  the 
altitudinal  direction.  Using  typical  resolutions  of  20  km  per  voxel  in  the  Y  direction,  the  smaller 
the  angle  of  elevation  of  the  nearest  point  on  the  orbit  from  a  ground  station,  the  greater  number 
of  longitudinal  planes  traversed  and  the  greater  the  cumulative  approximations  in  the  estimation 
of  the  discretized  altitude  of  features  as  shown  in  Figure  9b.  In  the  absence  of  corrective 
influences  from  the  sub-orbital  stations,  the  combination  of  these  two  factors  is  responsible  for 
the  upward  shift  of  the  band  in  the  reconstruction  of  Figure  8c. 

Figures  7c  and  8d  demonstrate  the  role  played  by  abruptness  of  sky  coverage,  due  to 
the  incomplete  recording  of  overpasses,  in  the  development  of  false  distribution  features.  The 
northern  half  of  the  coverage  map  is  affected  by  the  narrow  angular  span  of  coverages  from  the 
displaced  stations  3  and  4.  Relative  to  the  smooth  and  horizontally  homogeneous  reconstruction 
of  the  southern  band  in  Figure  8d,  the  northern  arm  appears  to  develop  variegated  regions 
resembling  enhancements  and  depletions  [Davies,  1990]  at  or  around  the  sites  of  abrupt 
coverage.  The  limitations  of  sky-coverage  is  compounded  by  the  non-coplanarity  of  the  stations 


Figure  9 
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in  this  region  which  adds  to  lack  of  information  in  the  relevant  projection  sets,  and  therefore,  in 
the  regional  density  distributions  of  the  reconstructed  image. 

Figures  8e  and  8f  correspond  to  reconstructions  of  data  sets  collected  using  the 
geometries  shown  in  Figures  7d  and  It.  The  first  imaging  system  consists  of  a  divergent  satellite 
orbit  and  receiver  chain  with  the  southernmost  station  as  the  only  sub-orbital  receiver.  In  the 
second  system,  the  geographic  disparity  between  the  satellite  orbit  and  the  chain  is  even  more 
pronounced  than  in  the  first  except  that  the  second  station  is  now  coplanar  with  the  satellite  orbit. 
Preliminary  logic  suggests  that  the  effect  of  the  greater  longitudinal  disparity  in  the  second 
instance  would  account  for  poorer  coverage  and  hence,  less  accuracy  in  the  reconstructions.  A 
comparison  of  the  southern  arms  of  Figures  8e  and  8f,  however,  prove  otherwise.  This  is 
explained  by  the  fact  that  central  location  of  the  coplanar  station  in  the  second  instance  enables 
the  corresponding  ray-path  fan  to  support  a  greater  area  both  in  the  latitudinal  and  altitudinal 
directions  whereas  the  support  of  the  coplanar  station  in  the  first  instance  is  greatly  weakened  by 
its  location  at  the  edge  of  the  chain.  These  two  geometries  are  frequently  encountered  in 
designing  CIT  systems  where  the  satellite  orbit  is  largely  over  vast  bodies  of  water  or  the  orbital 
plane  otherwise  inaccessible  for  receiver  deployment.  Under  such  circumstances  where  a  single 
coplanar  receiver  may  be  deployed,  better  results  are  obtained  by  positioning  it  towards  the 
center  of  the  chain  and  enabling  wider  sky-coverage  support  from  the  chain. 

5.  Conclusions 

The  preliminary  experimentation  with  coverage  and  image  resolution  for  ionospheric 
tomography  has  provided  valuable  insights  into  the  quality  and  quantity  of  information 
disseminated  in  the  course  of  image  reconstruction.  The  definition  of  coverage  used  in  these 
investigations  was  a  simple  and  practical  estimate  of  regional  availability  information  in  the 
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visible  portion  of  the  atmosphere.  It  is  anticipated  that  further  refinements  and  usage  of  this  tool 
will  include  other  tomographically  relevant  attributes  such  as  angular  distribution  of  ray-paths 
passing  through  a  voxel.  Resolution  analysis  has  demonstrated  a  simple  yet  elegant  procedure  of 
determining  optimal  computational  resource  allocation  to  meet  the  changing  requirements  of  this 
limited  angle,  limited  data  imaging  system.  The  development  of  some  commonly  encountered 
coverage  related  artifacts,  which  are  otherwise  difficult  to  identify  given  the  nature  of  the 
fundamental  image  features,  have  been  traced  to  specific  non-ideal  situations  created  in  the  sky- 
coverage  of  the  imaging  system.  It  is  anticipated  that  the  quantitative  basis  for  parametric 
estimation  of  sky-coverage  introduced  in  this  paper  will  encourage  further  investigation  and 
understanding  of  the  dissemination  of  information  during  the  data  acquisition  and  reconstruction 
process. 

The  procedures  of  coverage  and  resolution  analysis  presented  here  are  in  no  way 
specifically  tailored  to  suit  the  requirements  of  CIT.  They  represent  an  extremely  generalized 
approach  to  the  understanding  of  information  acquisition  and  distribution  in  geometrically 
restricted  tomographic  imaging  systems,  and  hence,  are  applicable  to  any  imaging  system  with 


similar  limitations. 
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Figure  4.  (a)  Typical  Variation  of  Percentage  Coverage  With  Image  Resolution  in  Two- 
Dimensions  (b)  Explanation  of  Coverage  Variation. 
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Figure  5.  Comparison  of  Coverage  Convergence  with  Resolution  for  (a)  MACE-93  and  (b) 
Australia-95  Campaigns  for  Z  =  4  of  7  Planes. 
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Figure  7.  (a)  CIT  Imaging  Geometry  and  Coverage  Map  for  System  1. 
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Figure  7.  (b)  CIT  Imaging  Geometry  and  Coverage  Map  for  System  2. 
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Figure  7.  (c)  CIT  Imaging  Geometry  and  Coverage  Map  for  System  3. 
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Figure  7.  (d)  CIT  Imaging  Geometry  and  Coverage  Map  for  System  4. 
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Figure  7.  (e)  CIT  Imaging  Geometry  and  Coverage  Map  for  System  5. 


Figure  8.  (a)  Original  Test  Image,  (b)-(f)  Reconstructions  for  CIT  Imaging  Systems  of  Figures 
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Figure  9.  Explanation  of  Altitudinal  Shift  for  Non-Coplanar  Stations,  (a)  Narrow  Angle  of 
Information  (b)  Low  Angle  of  Elevation. 


Coverage  Parameter 

Symbol  Type 

Definition 

Orbit  disparity 

A 

SOP 

Angular  divergence  of  satellite  orbit  from  center  longitude 

North  angular  fan 

CCx 

SOP 

Angular  coverage  to  the  north  from  station  x 

South  angular  fan 

A 

SOP 

Angular  coverage  to  the  south  from  station  x 

Fan  balance 

K 

SOP 

=  ( ax  -  Px )  /  (a,  +  P  J 

Chain  disparity 

K 

GSP 

Angular  divergence  of  chain  from  center  longitude 

Latitudinal  Spacing 

GSP 

Vector  of  receiver  spacing  in  north-south  direction 

Longitudinal  Spacing 

GSP 

Vector  of  receiver  spacing  in  east-west  direction 

Orbit  elevation 

GSP 

Elevation  of  nearest  point  on  orbit  w.r.t.  station  x 

Orbit  azimuth 

Xx 

GSP 

Azimuth  of  nearest  point  on  orbit  w.r.t.  station  x 

Directional  matrix 

A 

GSP 

Xl  »  X2  . . .  ^n>  Xn] 

Table  1. 


List  of  Coverage  Parameters. 
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#  Voxels 

MACE-93,  ICMT2-96 

Australia  -  95 

Caribbean  -  97 

[10, 60]  °N,  [15,65]  °N 

[-70,  0] 

[-10,  55]  “N 

°  Latitude  Earth-Dist 

°  Latitude 

Earth-Dist 

°  Latitude 

Earth-Dist 

/voxel  (km)  along 

/voxel 

(km)  along 

/voxel 

(km)  along 

X 

per  voxel 

X  per  voxel 

X  per  voxel 

10 

5.000 

555.0 

7.000 

777.0 

6.000 

666.0 

20 

2.500 

277.5 

3.500 

388.5 

3.000 

333.0 

30 

1.667 

185.0 

2.333 

258.9 

2.000 

222.0 

40 

1.250 

138.8 

1.750 

194.3 

1.500 

166.5 

50 

1.000 

111.0 

1.400 

155.4 

1.200 

133.2 

60 

0.833 

92.5 

1.160 

129.5 

1.000 

111.0 

70 

0.710 

79.3 

1.000 

111.0 

0.857 

95.1 

80 

0.625 

69.4 

0.875 

97.1 

0.750 

83.3 

90 

0.555 

61.7 

0.777 

86.3 

0.667 

73.9 

(a)  Latitudinal  Resolution. 

#  Voxels 

Altitudinal  Resolution  (km) 

[0,  1000]  km 

10 

100.0 

20 

50.0 

30 

33.3 

40 

25.0 

50 

20.0 

60 

16.7 

70 

14.3 

80 

12.5 

90 

11.1 

(b)  Altitudinal  Resolution. 

#of 

MACE-93 

Australia-95 

ICMT2-96  Caribbean-97 

Z 

[-105,-90]°E 

[140,  160]°E 

[-95,  -70]°E 

[-85,  -55]'’E 

plane 

°IZ 

Distance  (km) 

°/z 

Distance  (km) 

®  /  Z  Distance  (km)  ®  /  Z 

Distance  (km) 

1? 

3 

5.00 

278.3  548.1 

6.67 

742.5  253.9 

8.33  391.9 

895.7  10.00 

638.5  1096.3 

5 

3.00 

167.0  328.9 

4.00 

445.3  152.3 

5.00  235.2 

537.6  6.00 

383.1  657.8 

7 

2.14 

119.1  234.6 

2.86 

318.4  108.9 

3.57  167.9 

383.9  4.29 

273.9  470.3 

9 

1.67 

92.9  183.1 

2.22 

247.1  84.5 

2.78  130.8 

298.9  3.33 

212.6  365.1 

ll 

1.36 

75.7  149.1 

1.82 

202.6  69.3 

2.27  106.8 

244.1  2.27 

144.9  248.9 

t  Longitudinal  Resolutions  at  the  North  (N)  and  South  (S)  ends  of  the  view  volume  are  indicated. 

(c)  Longitudinal  Resolution. 

Table  2.  Chart  of  Geographic  Distances  and  Image  Resolutions  for  Geocentric  Coordinate 
System. 
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#  Z  plane 

MACE-93 

Australia 

ICMT2 

Caribbean 

1  of  3 

27-31 

0 

0-14 

2-4 

2  of  3 

66-74 

52-64 

55-90 

80-85 

3  of  3 

0 

0 

0 

0 

1  of  5 

0 

0 

0 

0 

2  of  5 

65-70 

C 

4-9 

48-56 

3  of  5 

30-35 

50-60 

72-78 

56-66 

4  of  5 

0 

0 

83-87 

3-6 

5  of  5 

0 

0 

0 

0 

1  of? 

3.5-6 

0 

0 

0 

2  of  7 

53-47 

0 

0-0.5 

10-17 

3  of  7 

63-52 

27-40 

11-16 

68-75 

4  of  7 

25-18 

44-54 

67-76 

67-74 

5  of  7 

0 

0 

86-88 

10-16 

6  of  7 

0 

0 

25-27 

0-1 

7  of  7 

0 

0 

0 

0 

1  of  9 

1. 4-2.8 

0 

0 

0 

2  of  9 

29-33 

0 

0 

2-4 

3  of  9 

52-61 

3-8 

2-6 

44-53 

4  of  9 

35-43 

28-38 

12-17 

48-58 

5  of  9 

11-16 

39-48 

52-62 

58-66 

6  of  9 

0 

0.5-2 

74-80 

14-20 

7  of  9 

0 

0 

59-65 

2-5 

8  of  9 

0 

0 

0 

0 

9  of  9 

0 

0 

0 

0 

1  of  11 

0.5-2 

0 

0 

0 

2  of  11 

18-23 

0 

0 

0 

3  of  11 

43-52 

0 

0-1 

14-21 

4  of  11 

42-53 

12-26 

6.5-9 

56-65 

5  of  11 

27-33 

23-33 

15-23 

31-41 

6  of  11 

7-11 

35-44 

48-58 

54-63 

7  of  11 

0 

2-6 

54-61 

16-24 

8  of  11 

0 

0 

78-85 

5-10 

9  of  11 

0 

0 

40-44 

0-2 

10  of  11 

0 

0 

0 

0 

1 1  of  1 1 

0 

0 

0 

0 

Table  3. 


Results  of  Percentage  Coverage  Versus  Resolution  Analysis  on  Four  Experimental 
err  Imaging  Systems. 
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Abstract 

This  paper  introduces  a  quantitative  feature  extraction  technique  for  the  evaluation  of  the 
content  and  the  quality  of  images  of  ionospheric  electron  distribution.  Parametric  alteration  of  the 
properties  of  image  histograms  of  electron  density  distributions  enable  the  separation  of  features 
such  as  enhancements,  depletions  and  small  scale  irregularities  of  such  images  from  the 
background.  This  technique  is  aimed  at  providing  comparison  of  plasma  images  obtained  from 
different  sources  and  for  the  utilization  of  the  image  information  content  for  various  applications 


in  real  time. 
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1.  Introduction 

Recent  advances  in  ionospheric  probing  and  imaging  techniques  and  in  the  development 
of  various  applications  which  utilize  images  of  ionospheric  electron  density  distribution  in  real¬ 
time  indicate  the  need  for  a  simple,  automatic,  unbiased  and  universally  applicable  numerical 
method  for  evaluating  the  contents  of  such  images.  Two  and  three  dimensional  maps  of 
ionospheric  electron  density  are  obtained  from  various  sources  ranging  from  extrapolated 
distribution  profiles  from  HF  ionosondes  [Davies,  1990;  Hausman  and  Nickisch,  1998],  to  true 
images  generated  by  ionospheric  tomography  using  VHF  and  UHF  signals  [Austen  et  al.,  1988; 
Fremouw,  1992;  Na  and  Lee,  \99A;  Raymund,  1994;  Fougere,  1995;  Na  and  Biswas,  1996;  Sutton 
and  Na,  1996;  Biswas  and  Na,  1998a;  Fehmers,  1998;  Kulinski,  1998]  and  even  more  accurate 
images  mapped  by  incoherent  scatter  radars  [Hunsucker,  1991].  The  geographical  extents  of  these 
images  can  vary  from  very  localized  areas  in  the  vicinity  of  quasi-vertical  ionosondes  to  several 
thousands  of  kilometers  below  a  satellite  orbit  or  along  the  look  directions  of  an  oblique 
backscatter  sounder.  The  resolution  and  reliability  of  the  actual  pixel  based  information  content  in 
these  electron  density  maps  vary  significantly  based  on  the  accuracy  of  measurement  and  on 
image  reconstruction  techniques.  In  addition,  the  information  content  of  a  given  image  is  defined 
by  the  applications  that  make  use  of  various  measurements  from  the  distribution  map. 

Current  literature  on  ionospheric  imaging  indicates  the  need  for  a  simple  and  unbiased 
numerical  method  for  evaluating  the  quality  of  reconstructed  images  which  is  suitable  for  diverse 
applications.  Images  obtained  by  computerized  ionospheric  tomography  (CIT)  have  traditionally 
been  subjected  to  visual  evaluation  and  analysis.  Comparison  of  CIT  reconstructions  with  images 
obtained  from  incoherent  scatter  radar  (ISR),  which  is  regarded  as  the  best  process  of  estimating 
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reconstruction  quality,  would  also  benefit  from  the  use  of  an  abstract  quantitative  analysis 
technique.  Often,  various  CIT  reconstructions  are  compared  to  each  other  for  algorithm 
evaluation.  This  situation  emphasizes  the  need  for  an  analysis  procedure  that  is  independent  of  the 
source  of  the  image  distribution  and  is  capable  of  comparing  image  information  at  a  more  abstract 
level  in  a  manner  similar  to  that  of  a  human  expert. 

Image  reconstruction  for  CIT  requires  the  simultaneous  participation  of  a  broad  spectrum 
of  disciplines  ranging  from  the  digital  domain  of  image  processing  to  that  of  plasma  physics 
involving  parametric  equations,  which,  for  all  practical  purposes,  can  be  considered  to  be 
functions  of  continuous  space  and  time  domains.  The  final  image  produced  by  a  generic  CIT 
algorithm  could  be  the  result  of  statistical  analysis,  iterative  algebraic  tomographic  reconstruction, 
basis  function  decomposition,  physical  constraint  based  optimization  or  various  combinations  of 
the  above.  In  addition,  the  scale  of  the  image  and  location  of  the  CIT  system  play  a  significant 
role  in  determining  contrast,  resolution  and  several  other  visible  attributes  of  features  within  the 
image.  Consequently,  the  list  of  image  differences  at  the  pixel  level  or  even  at  that  of  regional 
texture  is  often  too  exhaustive  for  enumeration  and  consideration  as  a  basis  for  quantitative  image 
analysis.  Various  real  time  applications  such  as  over-the-horizon  radar  (OTHR)  require  specific 
information  from  images  of  ionospheric  plasma  distribution  such  as  the  height  and  electron 
content  of  maximum  electron  density  [Biswas  and  Na,  1998b].  Pixel  based  comparative 
processing  throughout  the  image  plane  such  as  square  error  estimation  would  involve  unnecessary 
computation  when  all  that  is  required  is  the  mutual  comparison  of  a  set  of  one-dimensional 
profiles  from  a  large  number  of  competing  reconstructions.  Abstraction  of  an  image  using  feature 
extraction  and  classification  techniques  can  provide  robust  compaction  of  information  while 


5 


providing  a  scheme  for  application  dependent  quality  estimation.  It  must  be  emphasized  however, 
that  standard  image  comparison  schemes  such  as  mean  square  error  images  and  image  gradient 
differentials  carry  vital  information  regarding  image  quality  and  must  not  be  completely 
disregarded  in  the  analysis  process. 

The  aim  of  this  paper  is  to  introduce  methods  of  extracting  relevant  information  from 
competing  images  as  well  as  the  differentials  of  their  pixel  based  attributes  for  quantitative 
evaluation  of  quality.  The  next  section  discusses  the  major  properties  of  a  typical  image  of 
ionospheric  plasma  and  distribution  features  of  interest  to  major  applications  such  as  OTHR  or 
HF  communications  systems.  Section  three  introduces  the  processes  of  parametric  histogram 
distribution  (PHD)  and  subsequent  feature  extraction.  The  following  section  involves  a  discussion 
of  various  methods  of  image-plane  based  comparison  of  plasma  images  which  can  provide  direct 
measurement  of  several  significant  changes  during  real-time  monitoring  of  ionospheric  structures. 
Section  five  summarizes  the  need,  the  approach  and  the  significance  of  the  newly  developed 
image  content  and  comparison  technique  and  discusses  issues  associated  with  future 
developments  in  this  area  of  research. 

2.  Features  of  Interest  in  Ionospheric  Electron  Density  Images 

The  first  step  in  the  information  extraction  from  a  signal  for  abstract  representation 
involves  the  selection  of  features,  which  provide  robust,  yet  reliable  accounts  of  its  vital 
attributes.  Information  about  the  nature  of  commonly  encountered  ionospheric  features  was 
obtained  from  various  sources  ranging  from  direct  and  reliable  measurements  using  incoherent 
scatter  radar  (ISR)  to  images  generated  by  empirical  atmospheric  plasma  models  such  as  the 
International  Reference  Ionosphere  (IRI)  and  the  Parameterized  Ionospheric  Model  (PIM). 
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.  The  simplest  model  of  the  ionosphere  is  that  of  a  horizontally  stratified,  spherical  shell 
surrounding  the  earth  at  between  of  approximately  100-450  km.  The  altitudinal  variation  of 
density  can  be  approximated  to  an  asymmetric  Gaussian  or  an  alpha-Chapman  profile  [Davies, 
1990]  with  the  sharper  gradient  below  the  height  of  maximum  electron  density.  Observed  in  a 
two-dimensional  ground-altitude  plane,  the  section  of  this  shell  appears  to  be  a  horizontal  band  as 
shown  in  Figure  1.  Embedded  within  the  band  are  regional  enhancements  of  extremely  high 
electron  content  known  as  peaks  caused  by  plasma  convection  in  the  vicinity  of  the  magnetic 
equator  or  due  to  sudden  electromagnetic  disturbances  such  as  solar  storms  in  the  auroral  regions. 
The  size,  shape  and  density  of  these  peaks  vary  with  location,  time  and  ambient  heliogeomagnetic 
activity.  The  presence  of  such  embedded  peaks  contributes  to  significant  horizontal  gradients  in 
the  surrounding  band.  Isolated  patches  of  equally  high  density  but  significantly  smaller 
dimensions  occur  as  sporadic  E  and  F2  formations  that  drift  through  the  ionospheric  matrix  due  to 
gradients  in  plasma  density  between  the  sunlit  and  dark  portions  of  the  upper  atmosphere. 
Although  smaller  in  size  and  lasting  for  shorter  durations  than  the  peaks,  these  patches  are  of  vital 
importance  to  communications  systems  since  they  result  in  fading  and  radio  scintillation.  Thus  the 
band,  the  peak  and  the  patch  were  selected  as  the  three  most  important  features  of  an  image  of 
ionospheric  electron  density. 

Another  attribute  of  importance  in  an  ionospheric  image  is  the  profile  of  maximum 
electron  density.  This  can  be  viewed  as  a  surface  in  the  three  dimensional  spatial  co-ordinate 
system  of  latitude,  longitude  and  altitude  indicating  the  geographic  location  with  the  color  of  a 
point  on  the  surface  serving  as  the  fourth  dimension  indicating  the  maximum  electron  density.  For 
conventional  two-dimensional  CIT  images  with  the  altitude  versus  latitude  co-ordinate  system. 
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the  maximum  density  plane  reduces  to  a  curve  in  three  dimensions  with  electron  density  serving 
as  the  third  axis.  Examples  of  the  latter  convention  are  shown  in  Figure  2,  where  the  three 
composite  plots  represent  three  longitudinal  planes  of  three  IRI-90  images.  Each  composite  plot 
shows  the  variation  of  //max(F’2)  against  the  latitude  and  the  altitude  hmuxiFz)  in  the  three  images 
for  a  single  longitudinal  plane. 

In  addition,  pixel  based  image  differences  such  as  squared  errors  between  two  images  and 
their  normalized  versions  provide  important  information  regarding  the  recovery  of  absolute  and 
relative  electron  density  values  within  the  image  plane.  The  most  important  information  that  can 
be  extracted  from  square  error  plots  are  occurrences  of  isolated  pockets  of  high  error  such  as 
caused  by  relative  spatial  shift  of  features  between  two  reconstructions  or  a  reconstruction  and  an 
ISR  image.  The  indeterminate  nature  of  the  shape,  extent  and  magnitudes  of  these  regions  make 
it  difficult  to  classify  them  into  much  beyond  pockets  of  magnitude  and  shape  errors.  Comparison 
of  corresponding  pockets  in  the  mean  square  and  normalized  mean  square  error  images,  however 
yields  information  regarding  the  efficiency  of  value  recovery  in  the  region  versus  shape  recovery 
given  identical  conditions  of  coverage  and  image  reconstruction  procedure.  This  will  help 
distinguish  between  procedures  that  emphasize  value  recovery  in  localized  regions  from  those 
which  enforce  recovery  of  the  shape  of  plasma  distributions. 

Several  applications  such  as  ROTHR,  require  information  about  local  gradients  in 
ionospheric  images.  Differences,  in  both  the  magnitude  and  angular  components  of  gradient 
images  obtained  from  two  different  sources,  highlight  the  relative  susceptibility  of  the  inversion 
schemes  to  discontinuities  in  coverage.  As  with  mean  square  error  images,  the  principal 
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information  carried  by  gradient  differential  images  lie  in  the  existence  and  magnitude  of  high 
error  regions,  which  can  be  termed  as  slope  pockets. 

Features  chosen  for  image  representation  and  comparison  as  described  in  above  are 
enlisted  along  with  suitable  attributes  in  Table  1.  For  bands,  peaks  and  pockets,  their  considerable 
extent  with  respect  to  the  image  area  make  them  eligible  for  boundary  detection,  calculation  of 
geometric  centroids  and  value-centroids  or  valtroids,  and  the  length  and  orientation  of  the  major 
axis  in  the  Cartesian  co-ordinate  system.  For  small  scale  but  significant  structures  such  as 
patches,  the  major  attributes  for  classification  are  location  and  density.  The  second  column  of 
Table  1  refers  to  the  base  image  from  which  the  feature  is  extracted.  The  abbreviations  “ase”  and 
“nse”  stand  for  “absolute  squared  error”  and  “normalized  square  error”  images,  whereas 
“abgrme”  and  “abgrae”  stand  for  “absolute  gradient  magnitude  error”  and  “absolute  gradient 
angular  error”  images  respectively. 

3.  Feature  Extraction  Procedure  for  CIT  Images 

Feature  extraction  for  bands  and  peaks  was  performed  using  separate  intensity 
reassignment  procedures  for  the  histogram  of  a  given  ionospheric  image.  The  examples  of 
ionospheric  histograms  in  Figure  3,  show  the  presence  of  a  significant  background  mode  followed 
by  a  relatively  low  uniform  ledge  representing  the  foreground  with  most  of  the  electron  content. 
The  aim  of  the  histogram  reassignment  processes  is  to  separate  the  foreground  part  of  the 
histogram  into  two  separate  entities,  the  band  and  the  peaks  within  it.  The  image  processing 
technique  for  achieving  uniformity  in  histogram  profile  is  termed  equalization.  Therefore  the 
reverse  procedure  of  separating  a  uniform  histogram  into  one  or  more  distinctive  modes  can  be 
termed  unequalization.  Contrary  to  the  closed  form  definition  of  histogram  equalization  where  the 
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goal  is  to  achieve  a  uniform  histogram  density  for  all  intensities,  this  reverse  process  requires 
parametric  specifications  for  controlling  the  position  of  occurrence  and  extent  to  which  mode 
separation  is  permitted.  The  parametric  histogram  distribution  (PHD)  approach  developed  for 
ionospheric  images  involves  a  three-path  process  which  simultaneously  separates  the  band,  the 
peaks  and  the  patches  from  the  background  information. 

The  first  step  in  the  detection  of  bands  and  peaks  is  the  definition  of  an  intensity 
transformation  function.  In  the  case  of  histogram  processing  for  peak  detection,  it  was  necessary 
to  separate  the  original  histogram  into  two  distinct  modes,  that  of  the  peak(s)  and  that  of 
everything  of  lesser  intensity.  This  was  achieved  by  a  space  translated  single  term  polynomial  of 
the  form: 

v  =  (m- 0.5 +  £)■"+ 0.5 

Here  u  and  v  represent  the  normalized  axes  of  the  original  and  the  transformed  histogram 
intensities  respectively.  Figure  4a  shows  the  nature  of  the  above  curve  where,  for  all  practical 
purposes,  -0.5  <  8  <  0.5,  and  1  <  y  <  9.  The  parameter  y  represents  the  degree  of  contrast  enforced 
between  the  peaks  and  the  rest  of  the  image,  whereas  e  represents  the  threshold  intensity  at  which 
the  contrast  is  centered.  Figure  4b  shows  the  curve  representing  histogram  transformation  for 

isolation  of  the  band.  The  aim  is  to  separate  the  histogram  into  three  principal  modes  representing 

the  background,  the 

band  and  the  peaks  while  maintaining  a  uniform  intensity  within  each  mode.  The  equation  for  the 
curve  shown  in  Figure  4b  is: 

0.5(2m/  m  e  [0.00,0.25] 

0.5(2m-1/  m  6  [0.25,0.75] 

Q.5{lu-if  M  €[0.75,1.00] 


v  =  ^ 


(2) 
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The  curve  parameter  5  has  a  value  in  the  range  [1,  9].  Figure  5  shows  an  original  ionospheric 
image  and  the  result  following  histogram  redistribution. 

The  second  step  in  feature  extraction  is  intermodal  boundary  detection.  Both  the 
transformed  histograms  shown  in  Figure  5  have  irregular  profiles  with  several  subsidiary  peaks 
and  valleys  in  each  principal  mode.  An  established  natural  method  of  threshold  detection  is 
iterative  nodal  propagation  [Biswas,  1993]  where  all  the  discrete  intensities  of  the  histogram  act 
as  nodes  in  a  neural  network.  The  network  is  fired  by  initiating  a  nearest  neighbor  propagation 
and  feedback  sequence  simultaneously  from  each  peak.  The  iterative  process  results  in 
modification  of  the  histogram  profile  by  accumulating  the  contents  of  arbitrarily  located 
neighboring  modes  under  the  principal  modes  and  smoothening  out  the  profile.  This  process  is 
continued  for  both  the  peak  and  band  histograms  until  two  and  three  modes  remain  in  them 
respectively.  The  intermodal  valleys  denote  modal  thresholds.  Figure  6  shows  the  flowchart  for 
the  computation  of  threshold  and  an  example  of  progressive  modal  thesholding. 

The  third  and  final  stage  of  the  feature  extraction  process  is  the  detection  of  features 
within  the  boundary  images  and  calculating  vital  attributes  of  each.  Boundary  detection  is 
performed  using  an  8-way  directionality  map  [Gonzalez  and  Woods,  1993].  Detection  of  false 
features  is  minimized  by  only  permitting  those  with  greater  than  a  certain  number  of  boundary 
elements  to  be  registered,  based  on  the  geographical  dimensions  of  the  voxels.  Figure  7  shows  the 
boundary  images  of  the  peak  and  band  images  shown  in  Figure  5. 

Following  peak  and  band  extraction,  three  important  sets  of  attributes  are  calculated  for 
each.  They  are  the  positions  of  the  geometric  and  the  value  density  centroids  [Gx,  Gy]  and  [\4,  Vy] 
and  the  length  and  orientation  of  the  major  axis  [Ma,  M,,].  Given  a  boundary  defined  as  the 
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ordered  and  directed  set  of  points,  B  =  {(xi ,  yO,  (xi ,  ^2),  •••,  (x„  ,  y„)}  and  the  original  image  P  = 
{p(xi,yj)},  the  centroids  are  defined  as  follows: 


n  n 


[G,,G^.]  =  H  Zx,,  Zy,] 


/  =  1  i  =  1 


n 


[V,,Vy]  = - - - [  Z  XiP(x^,yi)+  Z>',pU,>3',)]  i  6  Feature 

/  =  ! 


(3) 

(4) 


The  length  and  orientation  of  the  major  axis  of  the  feature  are  defined  using  conventional 
Cartesian  coordinates,  where  given  the  end  points  [ax,  ay]  and  [bx,  by],  the  major  axis  vector  is 
defined  as: 


=  [^(a, -by  +  {a^ -b/  ,  90/7r*tan\(a^-b^)/ia^-b^))]  (5) 

Effectively,  these  three  vectors  serve  as  the  coordinates  of  a  six  dimensional  feature  space 
where  each  feature  Fj  is  represented  by  the  sextuplet:  Fj  =  [Gxj,  Gyj,  Vxj,  Vyj,  Maj,  M„j].  Thus  F  can 
be  defined  as  the  domain  of  the  primary  feature  space  of  the  contents  of  the  original  image.  Figure 
8  shows  the  set  of  feature  attributes  calculated  using  the  parametric  histogram  distribution 
technique  for  the  two  images  of  Figure  3. 


4.  Auxiliary  Differential  Measurements  for  Image  Comparison 

The  feature  extraction  process  can  be  extended  to  the  realm  of  differentials  such  as  root 
square  error  and  gradient  magnitude  and  gradient  angle  images.  As  summarized  in  Table  1,  the 
purpose  of  analyzing  the  root  mean  square  error  image  resulting  from  two  corresponding 
ionospheric  distributions,  is  to  identify  regions  and  magnitudes  of  differences  in  absolute  value. 
Feature  extraction  from  the  error  images  results  in  compaction  of  important  information  such  as 
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the  location,  extent  and  magnitude  of  differences  between  the  two  images.  The  features  obtained 
from  various  error  maps  such  as  those  between  an  ISR  image  and  CIT  reconstructions  using 
various  algorithms,  can  then  be  compared  on  a  quantitatively  abstract  basis  to  see  which 
algorithm  truly  succeeded  in  incurring  fewer  artifacts  and  achieving  good  reconstruction.  Figure  9 
shows  the  important  differential  images  associated  the  two  ionospheric  images  of  Figure  3. 

The  raw  square  difference  image  shows  the  actual  degree  of  mismatch  in  electron  density 
between  two  images  and  is  useful  for  a  direct  comparison  of  image  values.  This  .is  useful  when 
comparing  a  reconstruction  with  a  standard  such  as  ISR  image  or  even  two  different 
reconstructions  derived  from  the  same  TEC  data  using  different  parameters  or  algorithmic 
approaches.  A  finer  analysis  of  the  success  of  each  reconstruction  requires  comparison  of  the 
relative  accuracy  of  feature  recovery  in  each  image.  Instead  of  comparing  actual  electron  density 
values,  it  is  necessary  to  compress  the  dynamic  range  of  each  image  to  lie  within  the  normalized 
range  [0,  1]  in  order  to  be  able  to  detect  relative  recovery  of  features  within  the  image  plane.  The 
principal  use  of  the  normalized  square  difference  image  is  to  identify  the  strengths  and 
weaknesses  in  the  reconstruction  processes  associated  with  each  image  in  dealing  with  issues  of 
non  ideal  system  geometry  and  incomplete  sky  coverage.  When  one  of  the  images  is  a  standard, 
this  difference  measure  would  prove  highly  useful  in  differentiating  between  artifacts  and  valid 
image  features.  For  instance,  the  square  difference  image  of  Figure  9a  shows  a  substantial 
difference  in  the  magnitudes  of  electron  density  over  the  region  extending  between  the  50'**  and 
the  90'*’  voxel  horizontally.  This  prominence  is  not  evident  in  the  normalized  difference  image, 
which  indicates  that  although  there  is  a  disparity  between  actual  electron  densities  in  that  region, 
relative  to  the  maximum  electron  densities  in  either,  the  area  is  of  nearly  identical  density  in  both 
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images.  If  the  two  images  were  reconstructed  using  different  methods  and  the  maximum  electron 
density  were  somewhat  different  in  either  image.  Figure  9a  could  lead  one  to  conclude  that  they 
were  divergent  representations  of  the  actual  electron  density.  On  the  other  hand,  Figure  9b  would 
help  reconcile  the  above  argument  because  the  shapes  of  features  recovered  in  the  inversion 
process  would  match  each  other  even  if  magnitudes  do  not. 

Differences  in  the  gradients  of  two  fairly  similar  images  provide  information  regarding 
relative  displacements  of  heights  of  maximum  density,  hmaxiFi),  and  also  of  the  differences  in  the 
density,  iVmax(F’2),  itself.  Figure  9c  and  9d  show  the  magnitude  and  orientation  of  the  gradient 
difference  between  the  two  test  ionospheric  images  of  Figure  3.  According  to  the  magnitude 
image  of  the  difference  in  gradients,  there  is  approximately  a  difference  of  3  voxels  in  altitude  in 
the  hmaxiFi),  between  the  two  images  for  the  left  half  of  the  image  which  eventually  gets  resolved 
in  the  other  half.  Compared  to  an  expert  judging  the  characteristics  of  electron  density  profile 
such  as  maximum  height  and  gradient  of  the  profile  by  directly  looking  at  the  original  images,  the 
gradient  calculation  approach  provides  one  with  a  simple  and  elegant  method  for  qualitative  and 
quantitative  evaluation  of  layer  height  and  density  differences  between  the  images.  The  role 
played  by  the  orientation  image  of  the  gradient  difference  is  that  highlighting  the  area  of  maximal 
gradient  discrepancy  between  the  two  images.  Since  the  vertical  profile  of  an  ionospheric  image  is 
often  steeper  and  more  dynamic  than  the  horizontal  profile,  differences  in  orientation  are  likely  to 
be  expressed  as  horizontal  lines  passing  through  the  region  in  between  the  /imaxCFijs  of  the  two 
images,  that  is  at  the  altitude,  where  the  difference  between  the  skyward  slope  and  the  earthward 
slope  of  the  two  images  is  maximum.  Since  the  orientation  map  tends  to  have  a  single  such  line  of 
difference,  given  the  unimodality  of  the  vertical  profile,  this  height  can  be  used  to  diagnose  the 
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principal  region  of  discrepancy  in  the  more  complicated  magnitude  map.  In  addition  to  providing 
^max(^2)  information,  the  magnitude  map  is  also  useful  in  locating  discrepancies  between 
horizontal  inhomogeneities  such  as  a  trough  between  peaks.  An  example  of  this  is  also  seen  above 
the  14*''  and  the  IS***  horizontal  voxel,  delineating  the  difference  in  trough  locations  between  the 
two  peaks  as  two  voxels  apart. 

5.  Conclusions 

Image  processing  experiments  were  conducted  on  various  IRI-90  images  with 
approximately  the  same  feature  distributions  and  the  results  of  feature  extraction,  classifications 
and  auxiliary  differential  calculations.  In  the  absence  of  the  availability  of  other  CIT 
reconstruction  algorithms  besides  MIVIA,  the  comparison  of  different  reconstructions  from  the 
same  TEC  data  was  not  attempted.  This  was  also  not  undertaken  with  the  understanding  that  it  is 
an  entirely  independent  and  vast  area  of  research  in  itself. 

The  principal  advantage  of  the  parametric  histogram  distribution  based  feature  extraction 
process  to  ionospheric  image  analysis  lies  in  its  quantitative  and  objective  approach  to  defining 
the  intensity  ranges  for  image  features.  Utilization  of  normalized  histograms  ensures  that  the 
features  are  selected  independent  of  actual  electron  density  values  and  the  process  is  universally 
applicable  to  all  ionospheric  images  regardless  of  geographic  extent,  ambient  ionospheric  activity, 
method  of  image  generation.  In  addition,  certain  image  difference  measures  have  been  analyzed 
for  the  extraction  of  important  data  such  as  shape  recovery,  difference  in  layer  heights  and 
gradients. 
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Figure  3.  Sample  Ionospheric  Electron  Density  Distributions  from  IRI-90  and  Histograms  of  Number  of 
Occurrences  vs.  Electron  Density  Values. 
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Figure  5.  Peak  and  Band  Separation  Using  Parametric  Histogram  Distribution. 
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Figure  7.  Boundary  Detection  of  Peaks  and  Band,  (b)  Detection  of  Small  Scale  Disturbances. 
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Figure  8.  Results  of  Feature  Extraction. 
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Figure  9.  Difference  Images  of  Ionospheric  Distributions  Shown  in  Figure  3.  (a)  Absolute  Square  Difference  Image  (b) 
Normalized  Square  Difference  Image  (c)  Magnitude  of  Gradient  Error  Image  (d)  Orientation  of  Gradient  Error 
Image. 
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Feature 

Image 

Attributes 

Significance 

Band 

original 

boundary,  centroid,  valtroid,  maj.  ax.  length  &  slope 

band  existence 

Peak 

original 

boundary,  centroid,  valtroid,  maj.  ax.  length  &  slope 

peak  existence 

Patch 

original 

position,  electron  density 

patch  existence 

ASE  Pocket 

ase 

boundary,  centroid,  valtroid,  maj.  ax.  length  &  slope 

value  recovery 

NSE  Pocket 

nse 

boundary,  centroid,  valtroid,  maj.  ax.  length  &  slope 

shape  recovery 

MAG  Pocket 

absgrme 

boundary,  centroid,  valtroid,  maj.  ax.  length  &  slope 

grd.  mag.  recovery 

SLP  Pocket 

absgrae 

boundary,  centroid,  valtroid,  maj.  ax.  length  &  slope 

grd.  orient,  recovery 

Maxdensity 

original 

3D  (4D)  plot  of  Ne  vs  latitude,  alttitude  (&Iongitude) 

Fi  critical  frequency 

Table  1.  Features  Selected  for  Representing  Contents  of  Ionospheric  Images. 


